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Redshifted 21-cm Signals in the Dark Ages 
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ABSTRACT 

We have carried out semianalytic simulations to build redshifted 21-cm maps in the 
dark ages. An entropy-floor model is adopted for planting protogalaxics in simulated 
minihaloes. The model allocates gas quantities such as baryonic mass and temperature 
to every A^-body particle and extensively exploits the particle nature of the data 
in the subsequent analysis. We have found that the number density of simulated 
minihaloes in the early universe is well described by the Sheth & Tormen function 
and consequently the signal powers of sim ulated minihaloes are fa r greater than the 
Press & Schechter prediction presented by iFurlanetto fc Oh! (|2006l h Even though the 
matter power spectrum measured in the halo particles at z = 15 is about an order of 
magnitude smaller than the intergalactic medium (IGM), the 21-cm signal fluctuations 
of haloes are, to the contrary, one order of magnitude higher than the embedding 
adiabatic IGM on scales, fc < 10 /iMpc -1 . But their spectral shapes are almost same 
to each other. We have found that the adiabatic signal power on large scales lies 
between the linear predictions of the infinite spin-temperature model (T s ^> T cm b) 
and the model with the uniform spin temperature equal to background value (T s = 
T^ 9 ). Higher preheating temperature (or higher background entropy) makes the power 
spectrum of signals more flattened because the hotter IGM signals are more thermally 
broadened and minihalo fluctuations dominating on small scales are more severely 
suppressed by the higher background entropy. Therefore, this model-dependent power 
spectrum slope measured on the scale of 100 ^ fc ^ 1000 /iMpc^ 1 will enable us to 
easily determine a best- matching halo + IGM model in future observations. 

Key words: early Universe - cosmology: theory - large-scale structure of Universe 
- diffuse radiation - methods: A^-body simulation. 



1 INTRODUCTION 

Until recently, the Large-scale Structures (LSS) of the uni- 
verse and the Cosmic Microwave Background (CMB) radi- 
ations have been two principal research areas in the astro- 
physical cosmology. They are distinct from each other in 
their nonlinearity and observation sc ales. The LSS (jjarrettl 
12004 iGott ' et al 1 l2005l : |Peebles|[l980l ) is a complex nonlin- 
ear structure forming relatively re cently (compared to the 
CMB) while the CMB radiation fPenzias fc Wilson! 1 19651 : 
ISmoot et all 1 19921 : ISpergel et ail l2005r i comes through the 
last scattering surface when the universe is still in the linear 
regime just 0.4 million years after the Big Bang. The LSS 
is a cumulative result of nonlinear gravitational evolution 
over the entire age of the universe (t ag e ~ 13.7 billion years) 
and, on the other hand, the CMB glow is a transient event 
when lights and baryons are decoupled from each other as 
the universe cools down below T ~ 3000K (we simply ig- 



* E-mail:kjhan(9!cita.utoronto.ca 
f E-mail:pen@cita.utoronto.ca 



nore the subsequent interactions of CMB photons with cos- 
mic objects such as evolving gravitational potential and free 
electrons when they travel from the CMB photosphere to 
us). 

However, between these two epochs, there is another era 
(11 < 2; < 1000) called the dark ages that have recently been 
regarded as a big cos mological reservoir of rich information 
on the early universe jLoeb fc Zaldarriaga|[200l : |Penll2004 
IFurlanetto et all l2009al : Bowmar] 120091 ') and could help us 
pin down cosmic parameters more accurately ( Coorav et alj 
|200S| ; iMao et alJbOQSl ; IFurlanetto et alj|2009bh . This epoch 
is essential to the study of the cosmology since it can cover 
the shortcomings of the CMB and LSS studies. Among those 
shortcomings are the difficulty of CMB observation on galac- 
tic scales for a more accurate determination of the power 
spectral index of matter field and the difficulties arising dur- 
ing the interpretation of LSS observations due to the non- 
linear evolution and complex biasing effects. 

In the dark ages, the density fluctuation is still in lin- 
ear or quasi-linear regimes and the universe remains dark 
by large because there seldom exist strong photon-emitting 
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sources such as stars or active galactic neuclei (AGNs) . How- 
ever, there is an astrophysical observable, the neutral hydro- 
gen, which may imprint its existence on the blackbody spec- 
trum of cosmic background emission. The neutral hydrogen 
emits or absorbs a photon at 21-cm as a bound electron flip- 
flopps its spin direction at the IS state. As the basic physics 
of radiative transfer is well known, we can easily decode the 
observed spectrum to get information on line-of-sight dis- 
tribution of the hydrogen quantities such as temperature, 
projected surface density, and line-of-sight velocity. Obser- 
vations of the abundant neutral hydrogen in the dark ages 
would, consequently, provide us with a powerful window of 
opportunity for obtaining wealthy physics on the IGM and 
birth of protogalaxies at the dawn of the universe. 

The emission and absorption of the redshifted 21-cm 
line on the background CMB spectrum indicate the inter- 
actions between the photon and neutral hydrogen: the spin 
temperature of a neutral hydrogen is governed by the CMB 
temperature, color temperature of Lya p hoton, and the 
kinetic temperature of the hydrogen gas dPurcell fe Field! 
1 19561 : lFieldlll959l : iFurlanetto et aill2006l : Eiratall2006l i Even 
after the decoupling epoch, the baryonic gas temperature 
is tightly coupled to t he CMB temperature until z ~ 300 
l|Furlanetto et aill2006l ) after which the gas temperature be- 
gins to drop more rapidly (T g ~ (1 + z) 2 ) than the CMB 
temperature (T cm b ~ (1 + z)). Because the gas density is 
still sufficiently high, the spin temperature of gas is coupled 
to the gas temperature through atom-atom collisions until 
2 ~ 100. However, as the gas density drops, p g ~ (1+z) , the 
hydrogen density is not any more sufficient to hold the spin 
temperature. As a result, after z ~ 25 the CMB radiation 
has played a dominant role in driving the spin temperature. 

There are several ongoing and forthcoming projects for 
observing the redshifted signals with the state-of-the-art in- 
terferometry tec hniques. The Gi ant Metrewave Radio Tele- 
scoped fGMRT. IPen et ail 120081 ) has been set up in India 
with 30 parabolic dishes of 45m diameter spanning 25km in 
a Y-shaped configuration for the target wavelength from 50 
to 1, 420MHz. One of main observational targets is the emis- 
sion from neutral hydrogen in the protogalaxies or protoclus- 
ters between r edshifts 3 and 10 . The Murchison Widefield 
Arra\FlfMWA. lLidz et al.l[2008l ) is designed for the observa- 
tion of the amplitude and slope of redshifted 21-cm power 
spectrum on scales, k ~ 0.1 - 1 /iMpc -1 , especially for the 
reionization epoch. T he LOFAE0 (Low Frequency ARray, 
IZaroubi fc s"ilk1l2005h will use the arrays of dipole antennas 
built across the European countries for the observations of 
redshifted hydrogen signals in the epoc h of reionization. The 
PrimevAl Structure Telescope^ (RaST. IPeterson et alf2004h 
will consist of log-periodic antennas in China targeting the 
first luminous objects in the epoch of reionization. And the 
Square Kilometer Arra£| (SKA) will be built in the south- 
ern hemisphere to map three-dimensional distributions of 
neutral hydrogen in the dark ages and reionization eras. 

The purpose of this paper is two-folded. Firstly, we 
want to apply an entropy-floor model to semianalytically 
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simulating protogalaxies in minihaloes (llliev et al. 20021: 
Furlanetto fc Leobll20o3: iMartel et al.ll2003l: IShapiro et ail 



20061 ) from iV-bodv particles. To exploit the particle nature 



of the simulation data, we implement a new method to ac- 
curately measure the optical depth even in highly dense re- 
gions. Secondly, by applying the power spectrum analysis to 
the generated maps we want to fully assess the minihalo con- 
tribution to the diffuse backgrounds. However, it is a much 
challenging job to observe individual minihaloes at the low 
frequency of redshifted 21-cm photons even in the recent 
future due to the small angular size (#haio < a few arcsec- 
onds) and the weak signals (T& < a few tens mK). Most 
of the ongoing projects have angular resolutions one or two 
orders poorer than needed to detect individual haloes. But 
the signals of minihaloes over the diffuse background can 
be detectable by the current planned interferometers and, 
moreover, the rapid evolving radio astronomy will make it 
possible to push back the current resolution limit beyond 
the minihalo scales in the future. Therefore, it is worthwhile 
to study the minihalo signals in the cosmic context. 

For this analysis, we run cos mic JV-body simulations 
and apply the entropy-floor model JPenlll999l : lKaiserj|l99ll : 
IVoit et al.ll2005l : IOstriker et alj|2003 ) to build baryonic con- 
tents (or protogalaxies) in virialized dark minihaloes. Also 
the peculiar velocity and thermal broadening are included 
in the method to simulate the observed redshift distortions. 
The equation of absorption along the line of sight is fully 
solved even in halo regions where neutral hydrogen is so 
dense that the exact measurement of optical depth is much 
more important than anything else for the bright sources. 

We perform a dark-age benchmark test for redshifted 
signals of neutral hydrogen at z — 15 when the most of gas 
still remains neutral. In this paper, we do not include the 
effects of Lya and ionizing photons which are crucial for the 
study of the reionziation era. The onset of reionization epoch 
is still unclear because it is still beyond current observational 
barriers. However, we are able to get a clue from observations 
of quas ars and CMB. The pre sence of the Gunn-Peterson 
trough |Gunn fc Petersonlll965l ) in the quasar spectrum im- 
plies that there are abundant neutr al hydrogens in the IGM 
beyond z ~ 6 (jBecker et ail l200ll ) which implies that the 
reionization process is completed after this redshift. Also 
from the WMAP observations, the angular power spectrum 
of CMB consistently favors the extended ionization process 
jDunkelv et al.ll2009h that the reionization started at z ~ 11 
and finished at z = 7 ijSpergel et al.ll2~007l ). Therefore, the 
selection of z = 15 is adequate for the benchmark test for 
the dark ages. 

The content of the paper is as follows: the basic physics 
related to the 21-cm emission and absoption are given in Sec- 
tion [2] We introduce a new and robust method to generate 
the redshifted 21-cm signal maps from the JV-body simula- 
tion particles in Section [3] In Section [4] we describe how to 
measure the optical depth and how to achieve the doppler 
and thermal broadenings with the iV-body particles. Section 
[5] briefly describes the simulation and halo findings. Section 
[6] presents the resulting maps of various models and the ef- 
fect of the thermal broadening and Doppler shift on the im- 
age of minihaloes. Also the power spectrum analysis is given 
in the latter part of this section. We conclude with several 
arguements and remarks in Section [7] Appendix [A] is de- 
voted to provide a quick look at the differences of the linear 
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power spectrum between various methods. We extensively 
show how to determine the initial redshift of the simulation 
in Appendix [B] 

In this paper, we assume a concordance ACDM (cold 
dark matter) cosmology co nsistent with the WMAP 5-year 
data jKomatsu et al.l 20091 ') . We set the current CMB tem- 
perature to be T cm b(0) = 2.725, and the helium mass frac- 
tion, Y p = 0.24 (|Schramm fc Turnerlll998l '). We also assume 
that t here is no Lya photons so that the Wouthuysen-Field 
effect (|Baek et alJ(200Sl ; lHiratall2006f ) is fuUy neglected. The 
frequency of the 21-cm line in a rest frame is ^21 = 1.42 
GHz. We mix the use of terms, haloes and minihaloes, for 
the same sense. Also baryonic matter means the gas mixture 
of hydrogen and helium. Gas temperature generally implies 
gas kinetic temperature. 



2 21-CM EMISSION AND ABSORPTION 
LINES 

The spin temperature (T s ) is a weighted sum of three com- 
peting temperatures such as the C MB (T* cm b), color (T a 



and gas kinetic temperatures (T„) J Purcell fc Field! Il95rj : 



iFurlanetto et alj|200fj| : IShapiro et all 12006) 



T a = 



l+Va+Vc 



(1) 



where the color temperature is related to the Lya photons. 
The collisional weighting coefficient is scaled to the CMB 
contribution as 



Vc = 



i;Cio 



(2) 



where C10 is the collisional de-excitattion rate 
jPurcell fc Field! Il956h and A 10 is the Einstein spon- 
taneous emission coefficient. We set y a = assuming that 
there are no stars and AGNs which could emit the Lya 
photons and ionize neutral atoms in the medium. 

The lower panel of Figure [T] shows distributions of spin 
and kinetic temperatures of the IGM as a function of the lo- 
cal density contrast at 2 = 15. Here the gas temperature is 
assumed to follow the adiabatic (thick) or isothermal (thin) 
processes. The adiabatic spin temperature moves toward the 
kinetic temperature as the density contrast grows and the 
isothermal spin temperature is dropping to the cold gas tem- 
perature which is below the CMB temperature at this red- 
shift. It is interesting to note that the overdense region could 
be colder than the CMB, provided 8 ^ 21. The spin tem- 
perature is not tightly coupled to the kinetic temperature 
even at high-density contrasts because the absolute density 
is not sufficiently high at z = 15. The upper panel shows 
the distribution of the spin temperature as a function of the 
kinetic temperature for the gas density contrasts, 8 P = 0, 10, 
100, 1000, and 10000. In the mean field (S p = 0), the spin 
temperature is nearly invariant of the kinetic temperature 
up to T g — 10 4 K because the gas density is much low. If 
the density contrast rises up to 1000, the spin temperature 
is almost coupled to the kinetic temperature. Therefore, we 
expect that the spin temperature in the virialized haloes 
may strongly be coupled with the gas temperature because 
of the high density. On the other hand, in the mean field the 
spin temperature is nearly invariant of the gas temperature. 
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Figure 1. Lower. Kinetic (short- dashed) and spin (solid line) 
temperature distributions for the adabatic (thick) and isothermal 
(thin line) IGM models at 2 = 15. The CMB temperature is 
marked by the dotted curve. Also we overplot the product of the 
spin temperature and optical depth (tT s ) in a long-dashed line. 
This quantity does not depend on the temperature model of the 
medium. Upper. Distribution of spin temperature as a function 
of kinetic temperature for various density contrast at z = 15. 



One can calculate the brightness temperature over the 
background CMB temperature (T cm b) as 



AT„(u) 



7-0) 



Ts(z') — Tcmb (2') -t'O) 



oIt'(v), 



(3) 



where ATb = Tb — T cm b, z' is the redshift of the source, T„ 
is its spin temperature, and t'(v) is the optical depth to the 
source at the frequency, v. The increase of opt ical depth due 
to a gas element at z' may be computed as l|Shapiro et al.l 

Eooi), 



dr(v) = 



3AgAi 



32tt Ts(z') H(z') 



nm(z') , , . , 



(4) 



where we have used the relations, d\n(u') — dln(l + z') 
and nm = (1 — Xhii)p 9 (1 — Y p )/mu- The gas density is re- 
lated to background gas density as p g = (1 + 5 g )g b g 3 where 
g b g 9 = p c (z)Qf, and p c (z) is the critical density at z. At 
z = 15, the mean ionization fraction of IGM hydrogen is 
Xmi M = 1. 92 x 10~ 4 which is obtained from the RECFAST 
package Q (|Seager at al.lll999T ). Hereafter, we uniformly ap- 
ply this RECFAST value to all hydrogens even though they 
are located in halo regions. This assumption is correct be- 
cause the virialized temperature and the gas density of a 
minihalo (Mh J; 1O 7 /i _1 M0) is not as high as to ionize the 
hydrogen atom through a collisional process. The ionization 
time scale is nearly infinite. 
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The Gaussian velocity dispersion spreads the profile, 
(i>), of emission or absorption lines along the line of sight 



1 



(v -v) = 



(5) 



where the observed (one-dimensional) frequency dispersion 
<t{u') is related to the kinetic temperature of gas as <r{y') = 
(y 1 J 'c)(kT g /mH) 1 ' 2 . Here, mu is the hydrogen mass. The 
combination of the cosmic redshift and Doppler shift in the 
frequency is expressed as 



l'21 



(v/c) 



1 + (Vc) 2 



(6) 



where v is the peculiar velocity toward the observer. In 
the mean field of no peculiar velocity, we can get r(0) = 
(3\%/32n)(A 10 n m (z)/H(z))(T i ,/T s (z)) by assuming that 
4>(y) = S(y) where 5{v) is the delta function. In the limit 
of r <C 1, equation © can be approximated to 



34xm(l + 5) 



0.023 



0.15 1 + z 
n m h 2 16 



1/2 



>(*) 



mK, 



(7) 



where we applied the approximation that H(z) ~ 
HoQmQ+z) 1 ' 6 . In Figure[T] we show the distribution of tT s 
which is equal to the unredshifted brightness temperature in 
the limit of r <S 1. 



used because of its compact form beneficial for the analy- 
sis of the inner structure s of observed clusters (|Voit et alj 
l2005l : iMitcheiT et alJliooih . During the adiabatic cosmic ex- 
pansion, the background entropy (-Kigm) is fixed with time 
because T g oc (1 + z) 2 and p g oc (1 + z) 3 . However, the halo 
entropy in the SIS model is a rising function with radius, 
Ksis oc r< 4 / 3 \ and, as a result, the entropy below a critical 
radius may happen be less than Kigm- 

There have been many ob servational evidences 
jVoit et al.ll2005l ; iBalogh et afll2006h for the entropy floor 
in the inner region of cluster haloes and many researchers 
have pro posed various models for descri b ing the flat core 
entropy ilPenl 1 19991; lOh fc Haimanl 120031; IXu fc Wul 120031 ; 
iRovchowdhurv et all [20041 : lOstriker et alj 120051 ). In this 
work, we adopt the entropy-floor semi-isother mal sphere 
(EIS) model of the density distribution given by |Penl (|l999h 
who proposed a density profile as 



Pair) 



= (- 



3/2 



if r < R c , 
otherwise, 



where R c is the core size. 

The temperature outside the core is assumed to be 
isothermal and simply measured by equation l[8]l. And the 
core radius can be derived by equalising the entropy at the 
core boundary (r = R c ) to the background entropy, Kigm, 
as T v p g 2/Z (R C ) = T b3 (z)g b g s ~ 2/3 (z), where T b g a (z) is the 
background gas temperature at 2. If combined with equa- 
tions |[S} and @, this relation leads to 



3 MINIHALOES 

3.1 Isothermal Models 

In the singular isothermal sphere (SIS) model, the baryonic 
density of a virialized halo follows a simple power law as 
p g (r) — v 2 Q,b/iTTGfl m r 2 , where the circular velocity is de- 
fined as v 2 = GM V /R V . The real-space virial radius (R v ) is 
defined by the extent to which the mean density of the halo 
is vnsPc(z). At hig h redshift such as z — 1 5, it is sufficient 
to Set 1)178 — 178 (|Brvan fc Normanlll998l ). We can derive 
the isothermal kinetic temperature of the baryonic sphere 
using the energy relation, 



l( 



(S) 



2 

where fin is the mean molecular weight of a gas mixture. 
The mean molecular weight is computed by the helium mass 
fraction, Y p , and the ionization of hydrogen atom, Xhii- 

However, the SIS model breaks down in the central re- 
gion where the second law of thermodynamics might be vio- 
lated: the central entropy happens to be less than the back- 
ground IGM entropy. According to the second law of ther- 
modynamics, the entropy of a system should only increase 
and, therefore, haloes forming out of the IGM should have an 
entropy value equal to or larger than the IGM entropy. The 
astrophysical entropy density is defined by K = T g pg~ 2 ^ 3 , 
which is different from the classical definitiorQ but widely 



7 the standard definition of entropy is related to the astrophysical 
conventional form as s oc 

In #3/2 



3 / fe B T g (0)(l + 
2 V pHmnG 



3/4 



"178 



Up. 



1/4 



= 1.808 



l + z 



16 



3/4 



r 9 fc9 (Q) 

0.0214A' 



3/4 



1.219 

M v 



-3/4 



10 4 h-^Mr. 



0.258 

1/2 



-1/2 



V 178 



We note that the relative core size is anti-correlated with 
the halo mass indicating that smaller minihaloes are more 
strongly affected by the entropy constraint. 

Also we can derive another relation between the halo 
mass and the infall baryonic fraction. Here, the infall bary- 
onic fraction is defined by the mass ratio of infall gas to 
the entire gas which is initially located in the collapsed "La- 
grangian volume" of the dark matter. Since the total gas 
mass in a halo is computed with M g — 4rr J Q fl " p g (r)r 2 dr, 
the infall mass fraction of gas is simply given by 



where 



M v 



1 - 0.581 (f^) if R C <R V , 



A 



otherwise, 



A / N 2 25 Jty 

A ^ S 125 6 a 



1-erf I -V25^ 



12 In a 



+ (31 - 12 In a) 725 - 12 In ( 

375 



(11) 



If f g = 1, all baryonic matter settles down to the halo centre. 
And f g = means that there has been no baryonic collapse 
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Figure 2. Mass dependence of relative core size, R c /R v (bot- 
tom), and infall baryonic fraction, f g (top panel), at z = 15, 20, 
and 50. Vertical bars mark the Jean's mass at those epochs. 



and, consequently, the halo has no galaxy or no gas in it. 
The gas-temperature profile of a halo is simply measured as 

2/3 



T„(r) = 



T g (0) 



T v 



Pg( r ) 
Pc(0)«i, 



if r < R c , 



otherwise. 



Figure [2] shows the dependences of the relative core size 
(bottom) and infall baryonic fraction (top) on the halo mass. 
The core size is an increasing function of redshift at a fixed 
halo mass and, consequently, the infall fraction f g is decreas- 
ing with redshift. We may check whether th e entropy-floor 
mode l contains the Jean's mass condition (|Gnedin &i Huil 
1 19981 ) which is defined as the mininum mass of a spherical 
overdense region that can gravitationally collapse overcom- 
ing the resistant thermal pressure. In the figur e, we note that 
the J ean's mass scale (calculated by Eq. 23 of lShapiro et alj 
120061 ) is roughly corresponding to the mass of a halo whose 
fractional infall-gas mass is f g ~ 0.1 - 0.2 (or the corre- 
sponding relative core size is (R c /R v ) ~ 2 - 3). This indi- 
cates that, on the Jean's scale, 80 to 90% of total baryonic 
mass inside the halo Lagrangian volume may not infall to 
the halo centre. So it is reasonable to think that the entropy- 
floor model may inherently have the Jean's mass criterion. 
Therefore, we may simply skip the setting of the minimum 
halo mass usually applied to the semianalytic process to 
build protogalaxies in minihaloes. 



3.2 Spin & Brightness Temperature Profile 

We investigate the distributions of spin temperature in the 
SIS and EIS halo models. The spin temperature is derived by 
equation (JTJ after we measure the gas density and tempera- 
ture in each halo model. Figure [3] shows radial distributions 
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Figure 3. The spin temperature profile for virialized minihaloes 
in the EIS (thick) and SIS (thin) model at z = 15. For clarity, we 
add the model name in parenthesis for haloes of mass Mj, = 100 
and 10 3 /i -1 Mq. The thin solid horizontal line marks the CMB 
temperature at the epoch. 



of spin temperatures for various halo masses. The SIS halo 
(thin) has a flat spin distribution in the inner region. And 
haloes of mass below M s ~ 3x 10 3 h~ 1 M.Q have a spin tem- 
perature lower than the background CMB temperature. In 
the EIS model (thick) the characteristic mass scale where 
the spin temperature is same as the CMB temperature, de- 
creases down to M s ~ 10 2 Ii~ 1 Mq. From this figure, we 
know that T s is approaching T cm b at the outer boundary of 
a halo mainly due to the low gas density: T s is decoupled 
from T g and coupled to T cm b through the Compton scatter- 
ing. In the inner part of the SIS halo the spin temperature 
is saturated to the uniform gas temperature. However, the 
spin temperature in a core of the EIS halo keeps rising for 
Mh 1O 3 M0 simply because the gas temperature increases 
with radius in the core. 

Integrating equation Q over v with an assumption that 
4>(v) = 5(v) leads to the radial profile of the optical depth 



r(u,r) 



( § 



0.023 



l + z 
16 



3/2 



0.15 



-1/2 



(12) 



In Figurefjjwe show the dependence of optical depth both on 
the halo mass and on the radial distance. In the SIS model 
(thin) the optical depth increases to the centre but it de- 
clines with halo mass at a given re lative radius (or at the 
same baryonic density) as noted bv llliev et al.l (|2002u . Also 
it is interesting to note that beyond the radius of 0.2Rn% the 
slope of the optical depth for the halo of M = 10 2 /i _1 Mq, 
is steeper than more massive haloes. This is because the 
radial distribution of y c (kinetic contribution to the spin 
temperature in Eq. {T}) becomes steeper due to the lower 
virial temperature. Therefore, the spin temperature (also the 
brightness temperature) approaches the CMB temperature 
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Figure 4. Profiles of the optical depth in EIS (thick) and SIS 
(thin curves) halo models at z = 15. Halo masses are written 
besides the curves of the SIS model. 



more rapidly in this low-mass halo. For more massive SIS 
haloes, the spin temperature is more tightly coupled to the 
isothermal gas temperature, so the profiles are parallel to 
each other: the distribution in the log-log scale shifts hor- 
izontally for different isothermal temperatures. The radial 
distribution of the optical depth follows that of the isother- 
mal density as r oc (1 + 5 g (r)) oc r~ 2 . The optical depth in 
the core of an EIS halo (thick) is not so steep as in the SIS 
model because the increase of gas density toward the halo 
centre is partially offset by the increase of spin temperature 
(see Eq.[T21). 

Because haloes are ususally optically thick as shown 
above, it is important to apply a full expression of optical 
depth to the brightness temperature as, 



AT b (is) 



T s — Tcmb(z) 



-tO) 



(13) 



It is valuable to note that if r 2> 1, the brightness temper- 
ature is proportional to (T s — T cm b)- On the other hand, if 
t « 1 and T s » T cmb , we expect AT b x(l + S g ). Figured] 
shows the predicted radial profile of the brightness tempera- 
ture (Tb — ATb+r C mb(0)) for various halo masses. In the SIS 
model (thin curves), the temperature profile has two phases; 
the inner flat stage which is optically thick (t> 1), and the 
outer power- law stage which is less thick (r < a few). As 
expected, we observe that a massive halo with a high spin 
temperature shows a power-law brightness temperature pro- 
file in the outer region because the baryonic density is a most 
dominant factor in determining the brightness temperature. 
Also in the EIS model (thick), the slopes of the brightness 
temperature and the spin temperature in the core region are 
almost same to each other because of the nearly flat optical 
depth which virtually fixs the second term of the right-hand 
side of equation (JT3J. Minihaloes of M = 100 h~ 1 M Q in 
the SIS model may be observed as cold spots but those EIS 
counterparts can not be distinguished from the background 
CMB temperature at z — 0. 
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Figure 5. Profiles of brightness temperature of virialized mini- 
haloes in the EIS (thick) and SIS model (thin curves) at z = 15. 
The horizontal line guides for the background CMB temperature. 
The halo of M = 100 h~ 1 Mq in the EIS model, has a temperature 
profile nearly overlapped with the background CMB temperature. 



3.3 Halo Contribution to Diffuse Backgrounds 

Now, we investigate the contribution of minihaloes to the 
diffuse backgrounds in terms of the brightness tempera- 
ture and observed flux. The averaged brightness tempera- 
ture over the volume of a halo is defined in comoving space 
as 



(AT b (M v )) 



= 3 



4tt J*" AT b (M v ,r,e,4>)r 2 dr 
V(M) 

f AT(M v ,s)s 2 ds, 
Jo 



(14) 



where Vh(M) is the comoving volume of a halo of mass M, 
s = rjR%, and R% is the comoving virial radius. If the sky 
is uniformly illuminated with a brightness temperature, Tb, 
the observed flux is Tf,Afi a Ai/ bs where AQ. a is the antenna 
beam solid angle and Av a ^ s is the observation bandwith. 
Therefore, the total contributon of minihaloes of mass M 
to the diffuse background flux is a product of a single-halo 
contribution with the mean number density of haloes at the 
given mass as 

AT b (M)Au ohs An a = (AT b (M v )} An h Av h 

where Afih is the observed halo solid angle, Avh is the ef- 
fective halo size along the line of sight in frequency, <f>(M) 
(= dn(M) / d\og 10 M) is the number density of haloes of 
mass M, Az = Az(A^obs), and V is the survey comoving 
volume. Using equation (|15|l . one may easily derive 

AT b (M) = (AT b (M)) V h (M)$(M). (16) 

Here, we have used dV/dzdQ, — cd 2 /H(z) and AQh = A/d 2 , 
where c is the speed of light, d is the comoving distance to 
the halo, and A is the geometrical cross section of the halo 
in the comoving space. The product of last two terms is 
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Figure 6. Brightness temperature of individual minihalocs (bot- 
tom), mean brightness temperature (middle), and flux (top) ob- 
served with A9 a = 10' as a function of halo mass M at z = 15 
(thick), 17 (intermediate), and 20 (thin lines). The solid and dot- 
ted lines are for the EIS and SIS models, respectively. The dashed 
curve show the EIS distribution computed by applying the Eisen- 
stein & Hu power spectrum to the abundance of minihaloes at 
z = 15. 

the total volume fraction of minihaloes of mass M. This 
result meets the reasonable expectation that the halo con- 
tribution to the brightness temperature should be a sim- 
ple product of one-halo contribution with the comoving- 
volume fractions of the haloes of the same mass. One may 
compa re this equation with the one expresse d in the real- 
space ijlliev et all 120021 : IShapiro et ail boOfJ ). The overall 
halo contribution to the diffuse backgrounds is measured 
by ATI = J °° AT b (M)dlog 10 M. Also, the observed average 
flux from minihaloes can be measured by (|lliev et al.ll2002l ) 

<Wv(M) = T^ T AT 6 (M)An a , (17) 

where we set AQ, a = ir(A6 a /2) 2 and A9 a is the simplified 
antenna beam angle. 

Figure [6] shows the distributions of three temperature- 
related observables as a function of minihalo mass: the mean 
brightness temperature of a single minihalo (bottom), the 
mean brightness temperature in a unit solid angle (mid- 
dle), and the mean minihalo flux (top panels) received by 
an antenna of A6 a — 10' at z = 15, 17, and 20. Each 
solid and dotted curves are for EIS and SIS models, re- 
spectively. The contribution of minihaloes to the diffuse 
backgrounds peaks around M = 4 x 10 h~ Mq and the 
peak position appears to be invariant of redshift. In the SIS 
model, there are negative contributions from minihaloes of 
mass M < 3 X 10 3 h~ 1 M Q while the EIS haloes always 
make positive contribution to the diffuse backgrounds. The 
model dependences can be more easily seen on the lower- 
mass scales because lower-mass haloes have relatively bigger 
cores where the gas density and temperature more seriously 




z 



Figure 7. The minihalo emission. Solid lines are showing the 
background diffuse radiation from minihaloes predicted with the 
ST function while dotted lines are based on the PS function. 
Thick, intermediate, and thin lines are estimated with power in- 
dexes, n s = 0.96, 1, and 0.92, respectively. 



deviate from the SIS models. To determine $(M) in this cal- 
culation, we have applied the power spectrum of the CAMB 
Source to the halo mass function of Sheth & Tormen (1999; 
ST). Dashed curves show the resulting effect of the Eisen- 
stein & Hu (1998; hereafter EH) power spectrum on the 
distributions at z — 15. Appendix[X]specifies the differences 
between these tw o power spectrum estimations. 

As shown bv llliev et alj (|2002h . we also check the red- 
shift distribution of the brightness temperature and the cor- 
responding flux emitted from our minihaloes. In the upper 
panel of Figure [7] the brightness temperature of minihaloes 
increases with time showing some significant deviations 
among different spectral indexes of n s — 0.96 (thick), 1 (in- 
termediate), and 0.92 (thin curves). Also the slopes of bright- 
ness temperature predicted from the Press & Schechter func- 
tion (PS; dotted curves) are steeper than the Seth & Tor- 
men. This is because the PS function underestimates mas- 
sive minihalo populations at high redshifts while overesti- 
mates the number density of less massive minihaloes at lower 
redshifts. The minihalo flux over an antenna of A8 a = 10' 
shows a similar distributions. 



4 APPLICATIONS TO SIMULATIONS 

In this section, we introduce a new Lagrangian scheme for 
building brightness-temperature maps in the dark ages. We 
show how to assign hydrogen gas density and gas temper- 
ature to each A r -body particle. We, then, exploit the parti- 
cle nature of the data to compute the optical depth and to 
generate distortion maps adding the effects of the peculiar 
velocity and thermal broadening. 
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4.1 Adiabatic Contraction in the IGM 

In the IGM, we have measured densities at a given position 
using an adaptive smoothing kernel to enhance the spatial 
resolution. The local density is measured with the 30 near- 
est neighbors by setting the smoothing length (h s ) equal to 
half the distance to the 30'th nearist neighbor. Then, we es- 
timate density at the positions of the IGM particles using 
the smoothing kernel, 



f (l-|g 2 + |g 3 )/(^ 
W A {q)={ (2-q) 3 /(4irh 3 s ) 




for < q < 1, 
for 1< q ^ 2, 
and otherwise, 



where q = r/h s . Under the assumptions of the adiabatic 
contraction and no additional heating sources, we can mea- 
sure the kinetic temperature of the IGM gas using T g = 
(T g (z)} (1 + 5 S ) 7_1 where (T g (z)} is the kinetic temperature 
of mean backgrounds, 8 g is the gas density contrast to the 
mean background, and 7 = 5/3 for a monoatomic ideal gas. 



4.2 Brightness Temperature of Simulated 
Particles 

The mean differential brightness temperature is discretized 
according to the finite volume element of frequency range 
(y — Av/2 ^ v ^ v + Av/2) and cross section AS as, 



AT h {D) = — 



1 

Av 

N(9) 



v+Av/2 



AT h {v")dv" 



v-Au/2 



= E 



T s (zj) — T cm b(zi) - Ti (v) 



1 + Zi 



A Ti (v) 



(18) 



where N{i>) is the number of particles lying along the line 
of sight on the cross section AS. The contribution to the 
optical depth by a single particle is 



An(P) = 



v+Av/2 
-Ay/2 



^) dv" 

dv 



3\%A w T*n m (zi) 
327rT„(zi)H(zi) 



rv + Av/2 
J v- Av/2 



(y" - Vi)dv".(V$) 



where Vi is the redshifted frequency of the 21-cm line emit- 
ted from a particle, i, and nm(zi) is the mean density 
contribution from the particle and is measured by tim = 
m(HI)p(A5Ad) _1 where m(HI) p is the neutral hydrogen 
mass of the particle and Ad is the spatial depth correspond- 
ing to the frequency channel width, Av. The optical depth 
to the i'th particle is simply a sum of Arj(v) for interven- 
ing particles (1 ^ j < i) between the observer and the i'th 
particle: 



(20) 



dr(v) 

3<i 



This computation may benefit from sorting and queueing 
N(v) particles with the distance from the observer. The dis- 
cretizations applied in e quations (1 1811 and Q2 0H are valid if 
At; is sufficiently small ( |Mellema et al.ll200q ). From equa- 
tion (|12l) . we have found that At ~ 0.03 at z = 15 and 
the gas particle remains optically thin in most cases for 
10 ^ z ^ 100 in the dark ages. 




AS 



Figure 8. A schematic side view of the line-of-sight optical 
depth. An observer is located in the far left side of the figure 
so that we can apply the plan parallel approximation. AS is 
the two-dimensional cross section and Av is the frequency width 
(channel) of the survey element. The optical depth is growing as 
moving to right side of the figure while the redshifted frequency is 
decreasing. The curves show the Gaussian distribution of dr„/dv 
(tx <p(v)T^~ H~ 1 (z)) for each gas particle. 



The Doppler sh ift by the pecular velocity is given by 
|Shapiro et al.ll2006t ) 



(21) 



where /3 = v/c and v is the line-of-sight peculiar velocity 
toward the observer. The last term in the right-hand side of 
equation (1191) was reserved for the thermal broadening and 
of a functional form as 

™+Ai//2 

4>{v" — Vi)dv" = 



-Av/2 

erf 



Az//2 - 



^2oAv) 



erf 



- Ai^/2 - i 
V2cn(v) 



(22) 



where <t; (= (vi / c){k,BT g / mu) 1 / 2 ) is the one-dimensional 
(line-of-sight) velocity dispersion. 

A schematic side view in Figure [8] illustrates how to 
measure the brightness temperature in a discrete volume 
element of a depth Av and width AS. The observer is as- 
sumed to be located at a far left side of the figure so that the 
condition of the plane-parallel approximation can be satis- 
fied. In the bottom panel, the circles represent simulation 
particles and central box region has a surface area of AS 
and frequence range Av. The line profile broadened by the 
thermal temperature (here we do not include the Doppler 
shift caused by the pecular velocity) are described in the 
top panel. Particles which do not substantially contribute 
to dT v /dv (and consequently to Tb) between v and v + Av 
are marked by open circles. The area of the profile is anti- 
correlated to the spin temperature T s {zi). The total optical 
depth is obtained by integrating the Gaussian profiles over 
the frequency range. 

It is reasonable to assume that the IGM particles carry 
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an equal amount of baryonic mass as m g = m p (fli,/n m ). 
But this simple relation does not hold any more in the halo 
region because the baryonic matter is typically decoupled 
from the dark matter through the hydrodynamic processes. 
During the infall it is subject to another force, the gas force 
from neighboring baryonic matter due to high densities and 
temperatures in haloes. As seen in equation ([SJl, the distri- 
bution of baryon matt er is usually different from the NFW 
profile of dark matter (jNavarro et all 1 19971 ) and, therefore, 
we should assign a different amount of baryonic mass to 
halo member particles according to modelled distributions 
of baryonic matter. First, we measure the matter profile from 
member particles of the simulated halo with several radial 
bins. Second, we calculate the density ratio of the baryonic 
matter to the simulated matter for a given model of SIS or 
EIS. Then, we can compute the amount of baryonic mass to 
be assigned to halo particles for each bin. 




M (h-'M ) 



5 iV-BODY SIMULATIONS & HALO 
FINDINGS 

5.1 Simulations &z Halo Findings 

We have upgraded the GOTPM (|Dubinski et al.l2004l by in- 
corporating the CAMB Soured (|Lewis et al.ll200fj| ) for gen- 
crating initial power spectrum. This upgrade is important 
in this study because the length scales of interest are very 
small (k > 10 h.Mpc^ 1 ) that the effect of baryons on the 
matter power spectrum is significant. For comparison, dif- 
ferences among various power spectra provided by various 
method are discussed in Appendix We adopt a cosmo- 
logical model consistent with the WMAP 5-year cosmology 
and set the initial redshift of the simulation Zi — 300 (for 
the reason of this choice, see the Appendix |B|) . We com- 
pute the linear power spectrum of combined matter (CDM 
+ baryonic matter) at z — 15 and linearly scale back the 
power amplitude to z — 300. This is because not only does 
the amplitude of matter power spectrum shift with redshift 
but also the spectral shape changes with time even in the 
early universe. Even though we are using the pure iV-body 
simulation, we want to obtain a simulated power spectrum 
of the combined matter at z = 15. We call it the core sim- 
ulation and list several characteristics of the simulation in 
Tabled] 

The friend-of-friend (FoF) halo findings are applied to 
identify virialized minihaloes with the simulated particles 
at z = 15. For the linking length, We employ the usually 
adopted value, 0.2d mcan - It is interesting to check whether 
the FoF mass function at high redshifts is well described 
by the ST or PS functions even on this small scales. Fig- 
ure [9] presents the mass functions of simulated FoF haloes 
(filled boxes) and corresponding analytic functions such as 
the ST (solid) and PS (dashed). A set of upper two curves 
is obtained by integrating the power spectrum over a com- 
plete range of wavelength while the other set of two curves 
is obtained by integrating the power over the wavelength 
confined to the simulation box. The box-size effect on the 
halo number density is clearly seen in the figure. But this 
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Figure 9. Effect of finite box size on the population of simu- 
lated minihaloes at z = 15. Each pair of lines are the analytic 
mass functions of ST (solid) and PS (dashed line) obtained by 
integrating over a complete range of the power spectrum (upper 
two) and over a range confined to the simulation box (lower two 
lines). Filled boxes are the mass functions of the simulation. 



effect does not matter in the power spectrum analysis which 
is the main topic of this paper. 

The simulated mass function is in good agreement with 
the ST predictions. This agreement is slightly different from 
previ ous results of many Lagrangian and Eulerian simula- 
tions dWise fc Abelll2008l ; llliev et al.|[2006l ; iReed et al.ll2007l ; 
iLukic et alj|2007h in which authors argued that they have 
detected 50% underpopulations of haloes compared to the 
ST predictions at high redshifts. 



6 TEMPERATURE MAPS & POWER 
SPECTRUM 

6.1 Brightness Temperature Maps 

As described in previous sections, we have estimated the ki- 
netic temperature and local baryonic density at the positions 
of halo and IGM particles, and have allocated these hydrody- 
namic quantities to iV-body particles so we can treat them 
as gas particles. Using these pseudo-gas particles, we run 
several semianalytic simulations targeted for quantifying the 
effect of semianalytic parameters. Among these parameters 
are the switches to turn on or turn off the signal sources such 
as minihaloes and IGM, and switches to add the Doppler or 
thermal distortions on the map. Also we measure the IGM 
temperature by chosing one of the adiabatic and isothermal 
processes. Table [1] summarizes the semianalytic models we 
have used in this paper. The naming convetions are as fol- 
lows: we use upper cases H or B if the halo or IGM particles 
are included in generating the map, respectively. And the 
trailing lower script denotes the temperature model of the 
IGM. Sometimes, we use trailing marks as (t), (p), or (tp) to 
denote that those models include thermal broadening, the 
peculiar-velocity distortion, or both of them, respectively. 
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Table 1. Simulation parameters 



N p 




^box 


Nstep 


Zi Zf 


h 




n b 


b m p e 


512 3 


512 3 


0.512 


1188 


300 15 


0.719 


0.96 0.258 


0.044 


0.742 1.26 71.6 h- 1 M G 0.1 fi^kpc 



Cols. (1) Number of particles (2) Number of grids applied to measure the Zel'dovich displacements (3) Number of steps (4) 
Initial rcdshift (5) Final redshift (6) Hubble parameter (7) Spectral index of P(k) (8) Matter density at z = (9) Baryon 
density at z = (10) Dark energy density at z = (11) Bias factor (12) Particle mass (13) Gravitational force resolution 



Table 2. Scmianalytic models 



Name 


Sif 


mal 


Temperature model 


Density model 


Halo 


IGM 


Halo 


IGM 


Halo 


IGM 


Had 


yes 


no 


EIS^, 


adiabatic 


EIS^ 




TTSIS 

Had 


yes 


no 


sisj 


adiabatic 


sis" 




HB ad £ 


yes 


yes 


EIS 


adiabatic 


EIS 


u 4 


Bad 


no 


yes 




adiabatic 




W 4 


H 20 


yes 


no 


EIS 


20 K 


EIS 




Hioo 


yes 


no 


EIS 


100 K 


EIS 




Hiooo 


yes 


no 


EIS 


1000 K 


EIS 




B 20 


no 


yes 




20 K 






B100 


no 


yes 




100 K 




W4 


B1000 


no 


yes 




1000 K 




Wi 


HB 20 


yes 


yes 


EIS 


20 K 


EIS 


11 4 


HB 100 


yes 


yes 


EIS 


100 K 


EIS 


u 4 


HB1000 


yes 


yes 


EIS 


1000 K 


EIS 




(HB)iooo 


yes 


yes 


1000K^ 


1000 K 


EIS 


W4 



Cols. (1) Model name (2) Halo contribution to signal map (3) IGM 
contribution to signal map (4) Halo temperature model (5) IGM 
temperature model (6) Halo density model (7) IGM density model 

a EIS halo model for temperature 
^ EIS halo model for density 
c SIS halo model 
^ the reference model 

e uniform halo temperature fixed to 1000K 



We call HBad a reference model and most of the compar- 
isons are made to this model. 

Two-dimensional projected temperature maps are 
shown in Figure [TO] for the reference model (HB a d, 
top-left) and three is othermal models investigated by 
iFurlanetto fc Obi (|2006h : the EIS halo models with isother- 
mal IGM of T g = 20 K (HB20, top-right), T g = 100 K 
(HB100, bottom-left), and T g = 1000 K (HB1000, bottom- 
right). In the HBad model the average IGM is colder than 
the CMB by about 0.7 mK and, moreover, part of IGM 
surrounding the overdense fila mnetary structures i s colder 
than average IGM as noted by IShapiro et al.1 (|2006l ). If the 
isothermal IGM is at T g = 20K, the observed IGM temper- 
ature is substantially lower than the background CMB tem- 
perature compared to the adiabatic case (see Fig. 1). But the 
IGM brightness temperature becomes higher once its tem- 
perature is higher than CMB temperature as can be seen 
in the bottom-left panel. As the background IGM temper- 
ature is raised, the halo signal is getting weaker and haloes 
become less visible. At T g = 100K, most of field minihaloes 
disappear and only massive minihaloes in crowded regions 
survive the hot IGM. According to the entropy-floor model, 
the higher background entropy makes a halo have a bigger 




Figure 10. Temperature maps projected along the line of sight in 
a cubic box of a side length 512 h~ x kpc in four representative EIS 
models (clockwise from upper-left panel, HB a( j, HB20, HB100O1 
and HBioo). Top-left panel shows temperature fluctuations of 
haloes and IGM using adiabatic backgrounds. Other panels show 
the differential temperature maps obtained by fixing the temper- 
ature of IGM to T IGM = 20 K (HB20, top-right), T IGM = 100 
K (HB100, bottom-left), and T IGM = 1000 K (HB1000, bottom- 
right). Observed temperatures are measured by averaging bright- 
ness temperatures along the line of sight over < L < 512 h~ kpc 
which corresponds to 88.7402 ± 0.0154 MHz at z = 15. 



but less dense core. This explains the weaker minihalo sig- 



nals in the hotter IGM. At T a 



1000K, the diffuse IGM 



is the dominant source of the observed signals showing hot 
complex structures around dense regions. Most of the hot 
signals (AT b > 50 mK) in the model of T g = 1000K are 
mainly coming from the hot IGM gas. 

Now, we show the effects of Doppler shift (or peculiar- 
velocity distortion) and/or thermal broadening on the red- 
shifted 21-cm map in Figure fTTI The x axis of the figure is the 
line of sight so the observer is assumed to be far left side of 
the figure. The peculiar and thermal distortions make haloes 
spread along the line of sight and their effects are especially 
significant in halo regions. In the distorted map, the total 
observed flux in the whole simulation box is substantially in- 
creased. This is because the amount of absorption is reduced 
in the distorted field as the heavily-obscured emission source 
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Figure 11. Spectral distribution of brightness temperature for 
the HB a( j model ignoring the effects of both the thermal broad- 
ening and peculiar shift (upper-left), and considering only the pe- 
culiar shift (HB a d(p), upper-right), only the thermal broadening 
(HB ac i(t), lower-left), and both of them (HB a( j(tp), lower-right). 
The x axis of the figure is along the line of sight with a frequency 
range, 88.7402 ± 0.0154 MHz. 



appears to be dispersed into the less optically thick region in 
the frequency space. The Doppler shift makes the map nois- 
ier than the thermal broadening in halo regions. Compared 
to minihaloes, the IGM experiences less distortions due to 
the lower temperatures and smaller peculiar velocities. 



6.2 Effects of Approximation on Power Spectrum 

To justify our semianalytic approach, it is crucial to com- 
pare our results with the well-known analytic solutions or 
with numerical findings given in other papers. And one of 
the most powerful comparisons is using the power spectrum 
analysis. We have measured the three-dimensional power 
spectrum on th e signal maps and have c ompared them with 
those given by IFurlanetto fc Oh! (|2006h who measured the 
signal power spectrum of the minihaloes based on the PS 
function under the assumptions of T 3 3> T cm b and r <C 1. 
For proper comparisions, we take the same approximations 
to equations (|18[) and (|19[) . It is worth noting that the con- 
dition of T s 3> Tcmb also satisfies that t < 1. In this 
comparison, the switches of peculiar-velocity and thermal- 
broadening are turned o ff. To obtain the nonlin ear power 
spectrum of minihaloes, IFurlanetto fc Ohl (|2006f ) have re- 
garded the Fourier transform of halo density profile (NFW) 
as the one-halo term of the power spectrum and have as- 
signed to each halo a uniform spin temperature measured 
at the mean gas density (/ 9 Wi7s) of the halo. In their anal- 
ysis, they used the PS function to populate minihaloes in 
their analysis. In the top panel of Figure 1121 the dotted 
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Figure 12. (top): The halo signal powers at z = 10 (thick) and 
20 (thin) for the IGM temperature models: adiabatic, T g = 20K, 
T g = 100K, and T g = 1000K, from the top most curves of 
the same styles. The d otted curves show the predictions of the 
IFurlanetto fc Ohl feOOd) while the solid curves are the measure- 
ments in this study. Also the two dot-dashed curves are the pre- 
dictions measured using the linear matter power spectrum and 
the approximation, (T a 2> 0), at z = 10 and 20. (bottom): Effect 
of the uniform spin temperature on halo-only powers at z = 15. 
The dotted curves are power spectra obtained under the assump- 
tions of a uniform spin temperature of a minihalo. The spin tem- 
perature is evaluated at the mean overdensity (pTi = /gfi78; 
IFurlanetto fc Ohll2006l) while solid curves are obtained by mea- 
suring the spin temperature for each member particle. From the 
top most in each set of curves we show the simulated power spec- 
tra of H a d, H20, H100, and H1000 models, respectively. We plot 
the H ac i model distribution with a thick curve. The dot-dashed 
curve is the linear power prediction measured under the approx- 
imations, T a 3> T cm b and t < 1. We also plot the background 
IGM level in long-dashed curve. 



curves show the signa l power of minihaloes reported by 
IFurlanetto fc Ohl (|2006t ) at z = 10 (thick) and z = 20 (thin) 
while two sets of solid curves with the same line thickness 
for the same redshift are computed from the numerical sim- 
ulation in this study. To get simulated halo data at z — 10 
and 20, we have run an auxiliary simulation with the same 
setting as the core simulation except the initial linear matter 
power spectrum directly measured at z = 300. To each red- 
shift data we employ four models: from the topmost curve, 
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Figure 13. Three dimensional power spectra of the 21-cm signals 
at z = 15. The linear power spectrum for the infinite spin temper- 
ature (T B 2> T cm b) and r 1 is plotted with a dot-short-dashed 
curve. Also, for comparison, the power spectrum measured by 
setting the spin temperature fixed to the mean background value 
is shown in the dot-long-dashed curve. The simulated nonlinear 
power spectra with infinite spin temperature are shown in solid 
gray lines. The dotted curves are obtained by assuming the negli- 
gible optical depth (r « 1). The thick solid black curves show the 
power spectra computed without applying any approximation to 
the measurement. Two long-dashed curves show the halo (thick) 
and IGM (thin) power spectra for T s S> T cmb . We add the power 
spectra of the (HB) 10 qq for comparison in short dashed curves. 
The left set of four simulated power spectra are obtained from 
another simulation of a box size Lbox = 128 ft _1 Mpc. 

the adiabatic (H ad ), T g = 20K (H20), T g = 100K (Hioo), and 
T g — 1000K (Hiooo) models. As can be see n, there are signif- 
icant deviations between the two methods. iFurlanetto fc Qhl 
l)2006l ) underestimates the power by a few factors at z = 10 
and by an order of magnitude at z = 20 with a bigger dif- 
ference in a higher isothermal IGM model. However, this 
difference of the signal power agrees with the fact that the 
PS predicts less haloes than the ST on massive scales and 
more massive haloes have higher power amplitudes. This dis- 
crepancy becomes larger at a higher redshift or in a model of 
higher background entropy because halo signals come from 
more biased objects. The dot-dashed curve shows the linear 
prediction for the matter field of an infinite spin temper- 
ature. In bottom panel, four dotted lines show the power 
spectra of halo signals measured in our simulations at z = 15 
assuming the uniform spin temperature (from the topmost 
curve, H ad , H20, H' 100 , and H' 1000 ) while a set of solid lines 
shows the corresponding power obtained by computing spin 
temperatures for each simulation particles. 

Now, our interests are shifted to the effect of the widely 
used assumptions, r <g 1 and T 3 3> T C mb, on the signal 
power. Note that the condition, T 3 3> T C mb, sufficiently sat- 
isfies t <C 1 at z — 15. We simplify the situation by ig- 
noring the halo model and measure the density and tem- 
perature at the position of every particle by the W4 and 
adiabatic assumption, respectively. To cover a wider wave- 



length scale (fc < lOft/Mpc), we run another simulation in 
a bigger box of a side length, Lbox = 128 h~ Mpc, starting 
from Zi — 80. Figure [T3] gives the resulting power spectrum 
of 21-cm signals computed by turning on or off the approx- 
imation, r < 1, and in various temperature models. The 
two kinds of dot-dashed curves are obtained from the linear 
matter power spectrum but with different approximations: 
the dot-short-dashed curve is directly evaluated by adopt- 
ing the approximations that T 3 T cm b and r <g 1 while 
the dot-long-dashed curve is measured by fixing the spin 
temperature to the mean background value ((T s 6s (15))). So 
we may expect the true adiabatic power spectrum should 
lie within these two boundaries in the linear regime. How- 
ever, the nonlinear clustering distorts the power spectrum 
by increasing the small scale powers significantly. The solid 
gray curves (Pnonimcar(T 3 > T cmb )) are tne simulated non- 
linear power for the limits, T s 3> Tcmb, and t < 1. Here 
we are able to clearly observe the nonlinear clusterings on 
the small scales fc > 30 /iMpc -1 . Around fc = 20 frMpc -1 , 
the small drop of simulated nonlinear power is nothing but 
the cosmic variance of the simulation. We call it a baseline 
model if no approximation is made to the spin tempera- 
ture and optical depth. The thick solid black curve is the 
power spectrum in the baseline model. On smaller scales 
(fc ^ 10 /iMpc -1 ) the slope is much steeper than the linear 

(Plmear(T s > T cmb )) and simulated (Pnonlincar(T s > T cmb )) 

models of the infinite spin temperature. However, the power 
shapes are similar to each other on larger scales. The dotted 
curve is plotted to show how significantly the signal power 
deviate from the baseline model by the assumption of r <C 1. 
The assumption of negligible optical thickness results in a 
rise to the small-scale power at fc > 10 /iMpc -1 but leaves 
no effect on the larger-scale signals. This scale-dependent 
deviation in the power spectrum may come from the scale- 
dependent fluctutations of optical thickness. On larger scales 
(fc < 10 ftMpc -1 ), the density fluctuation does not develope 
so much that the optical depth is generally negligible while 
small-scale structures are well developed and may have sub- 
stantially higher optical thickness. 

It is valuable to note a change of the amplitude and 
slope of power spectrum when an assumption on the spin 
temperature is applied. For comparison, we overplot the 
power spectrum of (HB) 1000 model (thick- dashed). This 
model is added to check whether the gas temperature of 
T g = 1000K would be enough to make the spin tempera- 
ture sufficiently high to satisfy the infinity approximation 
(T s T cm b). The model predicts a much lower amplitude 
by a factor of about five than the P n oniinear(T s 3> T cm b) 
model but with the same shape. This amplitude difference 
may be caused by the insufficient gas density which finds it 
hard to efficiently pump up the spin temperature at z = 15. 
Even though the gas temperature reaches T g — 1000K, the 
spin temperature of a region of S p = (10) could only rise 
to T s ~ 50K (100K). From the matter power spectra of the 
halo and IGM, we measure the signal power spectra under 
the approximation of T s ^> T cm b. In the figure, the halo 
power (thick long-dashed) has a steep slope and crosses the 
IGM power (thin long-dashed) at fc ~ 6 x 10 2 /iMpc -1 . We 
do not draw the halo power in the bigger-box simulation 
because no virialized haloes (Mh 3.4 x 10 10 ft _1 M0) are 
identified due to a low mass resolution. 
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Figure 14. Halo power spectra in SIS and EIS models. The solid 
line shows the halo power of the entropy model (H ad ), and the 
dotted line is for the SIS haloes (H^ s ). For comparison, the power 
spectrum of adiabatic IGM (B at j) is shown in the dashed curve. 
The dot-dashed curve is the linear prediction for T s T cm t, and 

T < 1. 



6.3 Power Spectrum of haloes and IGM 

We want to highlight the effect of the entropy on the power 
spectrum of halo signals. The background IGM signals are 
excluded in order to isolated the role of entropy in the halo 
signals. Figure [T4l shows the signal power spectra of haloes in 
the two models (SIS in dotted and EIS in solid curves). The 
model difference is a function of the scale and increases with 
wavenumber, k. The EIS model shows a slightly higher am- 
plitude than the SIS model because of the higher core tem- 
perature and lower gas density which produces higher spin 
temperature and lower optical thickness in haloes, which 
leads to higher signal fluctuations (see Eq. I13p . The back- 
ground power of adiabatic IGM (dashed) has an amplitude 
an order of magnitude lower than all the available halo mod- 
els but has a similar power slope. 

It is important to note the characteristic scale below 
which halo signals begin to dominate the diffuse back- 
grounds. Figure [15] shows the power spectra of the halo 
(dashed), IGM (dotted), and both of them (halo+GIM solid) 
in three different halo models. In the HBiooo, the power spec- 
trum is dominated by the IGM and halo signals are com- 
pletely buried in it. Therefore, if the IGM was preheated to 
T g = 1000K, we can not observe the minihaloes in the red- 
shifted 21-cm observations. Meanwhile, the HBioo model has 
a dominant IGM power on larger scales (k < 50 feMpc -1 ). 
But, on the other hand, the halo signal is stronger on the 
smaller scale. In the HB a d, the halo power is much higher 
than the IGM showing the strongest power compared to 
isothermal models. Preheating before minihalo formation 
may be an inportant factor to expect whether minihaloes can 
be observed in the dark ages. A higher temperature of IGM 
more strongly suppresses the halo signals, leaving hotter and 
bigger cores of protogalaxies in minihaloes. The power am- 



Figure 15. The contrast of minihalo signals in the background 
IGM emission. The power spectra of redshifted 21-cm signals of 
minihaloes (dashed) and IGM (dotted) are shown against that of 
combined signals (solid curves) for HBiooO: an d HB a( j models. 
Each line is tagged with the applied model name. 



plitude of the isothermal IGM is an increasing function of 
the gas temperature and the power spectrum is steeper than 
the adiabatic IGM. Therefore, the slope of combined power 
spectrum depends on the IGM temperature at the epoch of 
minihalo formation. A steeper slope of the power prefers a 
lower IGM temperature and by observing the slope of the 
power at these scales we will know the preheating history 
of the IGM. The adiabatic IGM has a more steeper slope 
(riA ~ 1) while the IGM preheated to T g = 1000K has a 
slope of n A ~ 0.5 over 10 < k sC 1000 /iMpc" 1 . 



6.4 Effects of Thermal Broadening and Doppler 
Shifting 

As easily seen on the temperature map in Figure 1111 the 
thermal and Doppler broadening are significantly stretch- 
ing hot halo signals along the line of sight. Rich structures 
in the crowded regions of haloes are significantly smoothed 
out and a large fraction of small minihaloes in the mean 
fields are severly buried in the IGM due to the flattening. 
In Figure [TB] we show changes of the power spectrum from 
the reference model owing to the peculiar velocity (short- 
dashed) and thermal broadening (dotted). Th e peculiar ve - 
locity tends to enhance the large scale powers (|Kaiser||l987T ) 
but lowers the small scale powers as typically observed in 
the g alaxy redshift surveys (|Heavens et al.|[l998l ; IPark et"all 
Il994h . The power spectrum of biased objects distorted by 
the expone ntial or Gaussian velocity dispersi ons is well de- 
scribed by (|Park et al.ll 1994 ICole et al.lll995h 



P(k,n) =P R (k) 



2\2 



(1 + k 2 a^ 2 /2Y 



P(k, fi) = P R (k) (1 + /3/j 2 ) exp (-k W/2) 



(23) 



(24) 
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Figure 16. Power spectrum distortions caused by the thermal 
broadening and peculiar shift at z = 15 in the reference model. 
The thick solid curve is measured when we do not apply any of 
the distortions. The dotted curve is obtained if only the thermal 
broadening effect is considered and the effect of the Doppler shift 
is shown in the short dashed curve. The long dashed curve is 
the resulting power spectrum when two effects are simultaneously 
considered. Two gray lines are the predicted power spectrum. The 
dark gray and light gray are predicted from the Gaussian and 
exponential distribution models, respectively. 



respectively. Here P R (k) is the real-space power spectrum, 
(3 ~ f2°' 6 (,z)/&2i, &2i is the bias factor for 21-cm signals, a v 
is the velocity dispersion, and [i is the directional cosine of 
the wave vector along the line of sight. The former equation 
has been known for better de scription of the d istribution 
of simulated peculiar velocity (|Park et al.lll994l ) while the 
latter one is better for the Gaussian velocity dispersion like 
the thermal broadening. By simply scaling the large scale 
power we obtain 621 = 0.8 since fi(z) ~ 1 at the redshift of 
interest. Also the one-dimensional velocity dispersions of the 
simulated particles give that a v = 2.48(1 + z)km s -1 H -1 (z) 
for the peculiar velocity and a' v — 2.13(1 + z)km s -1 H -1 (,z) 
for the thermal broadening. Here, the redshift term is mul- 
tiplied to change the scale from the real space to comoving 
space. We integrated above equations over fj, and measured 
the averaged power spectrum. The thick gray curves in the 
figure show those predicted power spectra whose amplitudes 
are boosted up by I/&21 to match for the global amplitudes 
of the simulated power. On the scale of k ^ 300 /iMpc -1 
the amount of reduced power is much larger than expected 
for both distortions and this discrepency between the ex- 
pectation and measurment may come from the change of 
the optical depth as argued in the previous subsection. 

Now we study the effects of redshift distortions on the 
halo and IGM fields separately. In Figure 1171 the power 
spectra of the two distorted fields are shown with the same 
line types as shown in Figure 1161 To this study, we apply 
the thermal and peculiar distortions together. As seen in 
equations (|23|l and (|24|). there is nearly no global ampli- 
tude increase in the IGM power because the IGM has a 
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Figure 17. Same as Fig. I15l but including the distortion effects. 
Power spectra of H a( j(tp) and HB a( j(tp) are almost overlapped 
with each other. 



zero bias (b = 1) and £l(z) ~ 1. But the decrease of IGM 
power on the small scale depends on the temperature. In 
the HB a( j and HB100 models, the A(fc) have a peak around 
k — 1400 /iMpc -1 while HB1000 model has a power spectrum 
peak at k = 80 /iMpc -1 . We learn that haloes are a more 
dominating factor in the signal power than the background 
IGM on minihalo scales. Also in the HB100 model haloes 
are more powerful on the scale below k ~ 30 /iMpc -1 . And 
the most of the power spectrum in the HB1000 model comes 
from the IGM signals. The distorted IGM power explored 
in this paper is always smaller than the linear prediction 
estimated by using T 3 3> Tomb- The slope of the power spec- 
trum is considerably flatter than the non-distortion case: 
the HB1000 model shows a nearly flat power spectrum on 
100 < k < 1000 /iMpc -1 while HB ad model predicts a spec- 
tral index of tia — 0.5 on the same scale. 



7 CONCLUSIONS & DISCUSSIONS 

We have proposed a new semianalytic method to map the 
hydrogen distributions in the dark ages based on the La- 
grangian data of the N-body simulation. One of the most 
favourable features of the method is that it adopts a robust 
way of the optical depth measurement. Also the entropy- 
floor model is applied to properly describe the temperature 
and baryonic density in minihaloes. 

By analysing the power spectrum of the generated 
maps, we learn that haloes in the entropy-floor model dom- 
inate the adiabatic IGM in 21-cm signals over the entire 
scale (10 ^ k ^ 3 x 10 3 /iMpc -1 ) available in the simula- 
tions. However, the signal fluctuations of IGM with T g = 
1000K significantly overwhelm the minihalo signals over all 
scales. And the model with the IGM of a small temperature 
T g = 100K predicts that the halo power spectrum is larger 
than the IGM on the scale of k > 50(30) /iMpc -1 on the 
distorted (non-distorted) maps. Because the power spectra 
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of the IGM and halo have different slopes to each other, we 
may observationally determine the preheating temperature 
of the IGM from the power slope on the minihalo scales. 
Moreover, the thermal broadening and the peculiar velocity 
distortions make the slope more flattened. 

The adiabatic power spectrum of the matter field is on 
large scales five times larger than the linear prediction for 
a spin temperature fixed to the background value as T s = 
(T s 69 ) at z = 15. But on small scales (k > 10/iMpc' 1 ), the 
power spectrum rises more steeply with k than the linear 
model and becomes even larger than the upper-bound linear 
model of T s > T cmb at k ~ 600 /iMpc -1 . 

Test measurements of power spectrum with several ap- 
proximations frequently employed in the literature have 
shown substantial deviations around the true values. We 
have found that the worst case is the adoption of the Press & 
Schechter function for the minihalo number density at early 
universe. The Press & Schechter function underestimates the 
number of massive minihaloes compared to the simulation 
results which, on the other hand, is well described by the 
Sheth & Tormen function. Also the uniform spin temper- 
ature produces a considerably larger power than the true 
value. Either of the peculiar velocity and thermal broaden- 
ing seriously affects the power spectrum making the signal 
fluctuations globally boosted while significantly suppressing 
the small-scale fluctutions. Although this effect of the red- 
shift distortion on the power spectrum is similar to that of 
the galaxy surveys, we found that the analytic exponential 
or Gaussian distributions gives a poor fit to the simulated 
power spectrum on small scales. One possible explanation 
for this difference may be the decrease of the number of ab- 
sorbed photons in the distorted field. The heavily obscured 
regions are stronly subject to the distortions and the pho- 
tons emitted in the region can be less absorbed by the in- 
tervening hydrogen atoms because wavelengths of emitted 
photons are shifted to neighboring region where the optical 
depth would be lower. As a resut, this effect may suppress 
the small-scale signals less efficiently than expected. 

There are two assumptions employed in our numerical 
method. First, a gas particle is assumed to have a volume not 
being overlapped with the other particles. To measure the 
optical depth along the line of sight, simulation particles are 
sorted and queued in order of distance from the observer. 
Particles are considered to have their own exclusive box- 
shaped volumes of equal density and they are stacked along 
the line of sight with the same cross section. This approx- 
imation allows us to build a simple expression for optical 
depth as shown by equation (|20|l . If particles are allowed 
to have smooth density profile and they can overlap with 
each other, the equation of optical depth would be compli- 
cate and the implementation would be difficult. However, it 
would not seriously change the conclusions arrived in this 
study because the discreteness effect on the particle volume 
on the optical depth may be small and only confined to 
the resolution scale. Also the method does not consider the 
Wouthuysen-Field by the Lya photons emitted from stars or 
AGNs, which is very important in the reionization epochs. 
In the future extension of the method, we will cover this 
issue. 

Second, we assume the simulated haloes to be spheri- 
cal. In this semianalytic method, we have allocated baryonic 
mass and temperature to the JV-body particles. To generate 



protogalaxies in minihaloes, we subgrouped member parti- 
cles with multiple spherical shells of an equal centre. Because 
we already know the radial density profile of baryons in the 
EIS or SIS model, we can determine how much baryonic 
mass should be allocated to each bin particle. However, a 
problem happens when the distribution of simulation parti- 
cles of a halo is aspherical, which is common in most haloes. 
In this non-spherical case, the derived shape of the proto- 
galaxy built by our method is also following the shape of 
the simulated minihalo. So the spherical assumption in the 
EIS model is not valid any more. Therefore, care must be 
taken to the fact that the resulting baryonic distribution of 
minihaloes could be aspherical. 
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APPENDIX A: SIMULATION POWER 
SPECTRUM 

Comparing the linear power spectra proposed in the litera- 
ture is an interesting task and computing an exact power 
spectrum on small scales is crucial for the study of for- 
mation and evolution of haloes especially in the early uni- 
verse. For the CDM power spectrum Eisenstein & Hu (1998) 
proposed a simple fitting function that has been widely 
used for its simplici ty, fast speed, and easy implementation. 
Later, the CAMB (iLewis et alJlboOOf ) based on the CMB- 
FAST (|Seliak fc Zaldarriagalll996l '). provides a fast (but not 
so much fast as EH method) and error-controllable pack- 
age for the power spectrum. There is a variant of CAMB, 
the CAMB Source, which calculates more completely the 
small-scale transfer function considering the baryon sound 
speed. Figure IA1I shows the differences of the CAMB and 
CAMB Source from the Eisenstein & Hu's fitting function. 
In the left panel shown are the WMAP 5-year power spec- 
trum linearly extrapolated to z = 0. There seems no obvi- 
ous deviation between them below k ~ 100 ftMpc -1 . But 
if plotted in the differential power spectrum, we are able to 
notice their differences more easily. The differential power 
spectrum defined by SP(k) = (P(k) — PEri{k)) / PEii{k) are 
shown in the right panel of the figure. The solid curve shows 
the difference of CAMB Source from the EH and the dot- 
ted curve is for CAMB. In this figure, we note that the 
EH model estimates the matter powers within a few per- 
cent accuracy on the large scale (k > 10 2 fe/Mpc). And 
intriguingly, the CAMB slightly underestimates the small- 
scale power compared to the CAMB Source and can not 
properly recover the periodic ripples either on small scales. 
The difference of the power amplitude between the CAMB 
Source and Eisenstein & Hu is, on average, less than 4% 
with a maximum difference of 7% around k ~ 0.08 /iMpc -1 
for 10 -4 ^ k ^ 0.3 /iMpc -1 . But for scales of interest in 
this study (10 ^ k 3 x 10 3 /iMpc -1 ), the power spectrum 
difference between the CAMBs and the Eisenstein & Hu's 
becomes significant with high k value. This difference is sig- 
nificant when measuring the power spectrum of matter or 
redshifted 21-cm signals. So it is very important to use the 
CAMB Source for generating the initial power spectrum for 
the minihalo studies in the early universe. 



APPENDIX B: SETTING INITIAL REDSHIFTS 
& ITS EFFECT ON THE HALO MASS 
FUNCTIONS 

The GOTPM code adopts a first-order Lagrangian pertur- 
bation scheme to generate initial condit ions. It is faster an d 
simpler than the second-order scheme (|Crocce et al.l [20061 ). 
But the initial conditions should be generated with much 
care to satisfy assumptions adopt ed in the s c heme . Like 
the definition in the appendix of iKim et all (|2008l ), the 
Zel'dovich redshift (zk) of a particle is the redshift when 
the displacement of a particle from its Lagrangian point is 
equal to the mesh spacing either in x, y, or z direction. 
Prior to setting up the initial conditions, particles are lo- 
cated at the mesh points and linear velocities are assigned 
to them according to the Zel'dovich approximations from 
which the initial displacements are calculated. Therefore, if 
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Figure Al. Matter power spectrum at z = 0. (Left panel): Power spectra of the Eisenstein & Hu (dashed), CAMB (dotted), and CAMB 
Source (solid) are shown based on the WMAP 5-year cosmology. (Right panel): Deviations from the Eisenstein & Hu are shown for the 
CAMB (dotted) and CAMB Source (solid curves). All the power spectra are generated for the WMAP 5-year cosmology. 




the displacements are greater than the mesh size, they can 
not reflect the mesh-size fluctuations of the Zel'dovich poten- 
tials in the first-order scheme. So the initial starting redshift 
should properly be less than the Zel'dovich redshifts. 

There is an alternative to this definition of the 



Zel'dovich redshift: the redshift when the relative separation 
between neighboring particles is zero in one of the three di- 
mensions. This type of the Zel'dovich redshift would be the 
proper one if one wants to avoid the situation when two adja- 
cent particles overshoot each other in the initial conditions. 
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This constraint on the starting redshift would be more strin- 
gent than the original one. We call this the particle-particle 
Zel'dovich redshift (z pp ) and the former one is renamed the 
particle-mesh Zel'dovich redshift (z pm ). In Figure IB1I we 
show the distributions of z pm (left) and z pp (right) for sim- 
ulation resolutions expressed in terms of the mean particle 
separation (d mC an). In the simulations we use 256 3 particles 
and the same number of mesh to measure the initial dis- 
placement from the generated Zel'dovich potential. There- 
fore, one should note the lack of large-scale power in these 
"gauge" tests and the underestimation of the Zel'dovich red- 
shifts are expected especially for the small d mea n cases. For 
the z pm distribution, all distributions show an almost same 
shape with a single peak around which a power-law rising 
and a sharp cut off with z pm are shown. And the z pp distri- 
bution has double peaks with drops at higher redshifts than 
z pp . 

One of the easiest and simplest ways to justify the se- 
lection of the initial redshift is to measure the abundance of 
simulated haloes at later epochs. Simulations starting with 
redshifts equal to or higher than a certain critical epoch 
should show the same halo mass functions, and simulations 
with lower starting redshifts may have deviations from the 
true distribution. Figure IB2I emphasizes the importance of 
z pm in determination of the initial redshift. For simulations 
with cubic boxes of side length, Lbox = 0.512 h~ Mpc and 
512 3 -size mesh, we select three chracteristic initial redshifts 
for comparison. First, Zt = 500 and z% = 300 are chosen 
from distributions of z pp and z pm , respectively. And then, 
Zi — 100 is added to contrast the lower starting redshift 
against above higher values. As can be seen in Figure iBll 
most of initial particles at z% = 100 are shifted larger than 
the mesh spacing. Using these initial settings, we run the 
simulations down to z = 15. The numbers of time steps 
are determined to satisfy that the maximum displacement 
of particles in a step be less than O.ldmcan which is also set 
equal to the force resolution. The total time steps from Zi to 
z = 15 are set nearly same (ranging from 1,063 for 2, = 100 
to 1,214 for Zi — 500 simulation) in the three simulations. 
It is interesting to note that at z = 18 and 20 the simula- 
tion with Zi = 100 has underpopulations of haloes (~ 50%) 
compared to other simulations while the abundance differ- 
ences narrows to a few percent level at z = 15. There is no 
obvious difference between the z% = 300 and Zi = 500 sim- 
ulations so we conclude that it would be better to use the 
2 pm distribution for the starting redshift of simulations. 
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Figure B2. FoF halo mass functions at z = 15 (left), z = 18 (center), and z = 20 (right). Open boxes show halo abundances from the 
simulation with a starting redshift, Zi = 100 while open circles and filled circles represent the mass functions for Zi = 300 and Zi = 500, 
respectively. The solid curve follows the ST function and the dashed curve marks the PS function computed with the power spectrum 
confined to the box of a side length Lbox = 0.512 ft~ 1 Mpc. 



